germination

library(rethinking)

Loading required package: rstan
Loading required package: StanHeaders
Loading required package: ggplot2
rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Loading required package: parallel
rethinking (Version 2.11)

Attaching package: ‘rethinking’

The following object is masked from ‘package:stats’:

    rstudent

library(brms)

Loading required package: Rcpp
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Loading 'brms' package (version 2.13.5). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: ‘brms’

The following objects are masked from ‘package:rethinking’:

    LOO, stancode, WAIC

The following object is masked from ‘package:rstan’:

    loo

The following object is masked from ‘package:stats’:

    ar

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble  3.0.3     ✓ dplyr   1.0.0
✓ tidyr   1.1.0     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
✓ purrr   0.3.4     
── Conflicts ──────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x tidyr::extract() masks rstan::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x purrr::map()     masks rethinking::map()

germ <- read_csv("../Dimensions/hydrothermal-all-species/data/light_round1_tall.csv") %>% filter(wps==0) %>%

Warning message:
Unknown or uninitialised column: `germ`. 

  select(pops, temps, total_seeds, germ, day, cumulative_germ)

Parsed with column specification:
cols(
  pops = col_character(),
  temps = col_double(),
  wps = col_double(),
  date = col_character(),
  total_seeds = col_double(),
  germ = col_double(),
  start_date = col_date(format = ""),
  census_date = col_date(format = ""),
  day = col_double(),
  cumulative_germ = col_double(),
  cumulative_prop_germ = col_double()
)

germ

what if need one event per row:

one_per_row <- function(df) {
  total_seed <- max(df$total_seeds, sum(df$germ))
  newdata <- tibble(id=1:total_seed, germ=0, day=max(df$day))
  df <- df %>% filter(germ>0)
  count <- 1
  if (nrow(df) > 0) {
    for (i in 1:nrow(df)) { # we look at each row of the df where germination occured
      for (j in 1:df$germ[i]) { # now update the newdata to reflect the germiantion of each seed
        newdata$germ[count] <- 1
        newdata$day[count]=df$day[i]
        count <- count+1 # count keeps track of which individual we are at in the new data
      } # for j
    } # for i
  } # if 
  return(newdata)
}

germone <- germ %>% group_by(pops, temps) %>%
  select(-cumulative_germ) %>% # not needed in this encoding (I think...in any case would need to be recalculated)
  nest() %>%
  mutate(newdata=map(data, one_per_row)) %>%
  select(-data) %>%
  unnest(newdata)

germone

##1

Is effect of temp ~linear on cum germ?

germ %>% filter(pops=="STDI", day==28) %>%

Warning message:
Unknown or uninitialised column: `germ`. 

  ggplot(aes(x=temps,y=cumulative_germ)) +
  geom_col()

no, so we can’t treat temp as a continuous linear predictor

germ.stdi <- germone %>% filter(pops=="STDI") %>% select(-pops)

Adding missing grouping variables: `pops`

germ.stdi

Try a censored lambda model

d <- list(germ=germ.stdi$germ, temps=as.numeric(as.factor(germ.stdi$temps)), day=germ.stdi$day)

Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
3: Unknown or uninitialised column: `germ`. 

m1.1 <- ulam(
  alist(
    day | germ==1 ~ exponential( lambda),
    day | germ==0 ~ custom(exponential_lccdf( !Y | lambda)),
    lambda <- 1.0 / mu,
    log(mu) <- a[temps],
    a[temps] ~ normal(0,1)),
  data=d,
  chains=4,
  cores = 4
)

starting worker pid=54244 on localhost:11220 at 11:44:48.270
starting worker pid=54258 on localhost:11220 at 11:44:48.609
starting worker pid=54273 on localhost:11220 at 11:44:49.300
starting worker pid=54288 on localhost:11220 at 11:44:49.607

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000109 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 7.5e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 7.4e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.74 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.294589 seconds (Warm-up)
Chain 1:                0.230801 seconds (Sampling)
Chain 1:                0.52539 seconds (Total)
Chain 1: 
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 6.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.280862 seconds (Warm-up)
Chain 2:                0.212282 seconds (Sampling)
Chain 2:                0.493144 seconds (Total)
Chain 2: 
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.310661 seconds (Warm-up)
Chain 3:                0.219301 seconds (Sampling)
Chain 3:                0.529962 seconds (Total)
Chain 3: 
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.323819 seconds (Warm-up)
Chain 4:                0.215416 seconds (Sampling)
Chain 4:                0.539235 seconds (Total)
Chain 4: 

precis(m1.1, depth = 2)

Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 

posterior

preddata <- expand_grid(temps=1:8, realtemps=sort(unique(germ.stdi$temps)), day=1:28)

Warning message:
Unknown or uninitialised column: `germ`. 

pred <- link(m1.1, data = preddata)
str(pred)

List of 2
 $ lambda: num [1:2000, 1:1792] 0.00439 0.00507 0.00342 0.00459 0.00509 ...
 $ mu    : num [1:2000, 1:1792] 228 197 293 218 197 ...

mu is average days to germ, lambda is rate parameter. neither change over time, of course, so having days doesn’t really make sense.

preddata$mu <- apply(pred$mu,2,mean)

Warning message:
Unknown or uninitialised column: `germ`. 

preddata$low <- apply(pred$mu,2,HPDI)[1,]
preddata$high <- apply(pred$mu, 2, HPDI)[2,]

Single temp. Don’t need day posterior

preddata <- expand_grid(temps=3)

Warning message:
Unknown or uninitialised column: `germ`. 

pred <- link(m1.1, data = preddata)
str(pred)

List of 2
 $ lambda: num [1:2000, 1] 0.0176 0.0212 0.0178 0.0241 0.0221 ...
 $ mu    : num [1:2000, 1] 56.9 47.3 56 41.5 45.3 ...

how to convert to probs? use pexp.

predprobs <- pexp(1:28,rate=pred$lambda[1])

Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 

plot(x=1:28,y=predprobs, type="l")

even thought it isn’t using day, including day in the prediction data frame will help me keep data in the correct format. Maybe.

preddata <- expand_grid(temps=1:8, day=1:28)

Warning message:
Unknown or uninitialised column: `germ`. 

pred <- link(m1.1, data = preddata)
str(pred)

List of 2
 $ lambda: num [1:2000, 1:224] 0.00439 0.00507 0.00342 0.00459 0.00509 ...
 $ mu    : num [1:2000, 1:224] 228 197 293 218 197 ...

predresults <- preddata %>%
  mutate(lambda=as.list(as.data.frame(pred$lambda)))
predresults

predresults <- predresults %>%

Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 

  mutate(probs=map2(day, lambda, ~ pexp(.x, .y)),
         mu=map_dbl(probs, mean),
         lower=map_dbl(probs, ~ HPDI(.)[1] %>% unlist()),
         upper=map_dbl(probs, ~ HPDI(.)[2]) %>% unlist())
predresults

predresults %>% select(-lambda, -probs) %>%

Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
3: Unknown or uninitialised column: `germ`. 
4: Unknown or uninitialised column: `germ`. 
5: Unknown or uninitialised column: `germ`. 
6: Unknown or uninitialised column: `germ`. 
7: Unknown or uninitialised column: `germ`. 
8: Unknown or uninitialised column: `germ`. 
9: Unknown or uninitialised column: `germ`. 

  mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
  ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
  geom_line() 

Add realdata:

stdi.plot <- germ %>% filter(pops=="STDI") %>% 

Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
3: Unknown or uninitialised column: `germ`. 

  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop.germ=cumulative_germ/total_seeds)

predresults %>% select(-lambda, -probs) %>%
  mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
  ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
  geom_line() +
  geom_point(aes(y=prop.germ), data=stdi.plot)

Poor!